Towards the phase diagram of dense two-color matter 
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We study two-color QCD with two flavors of Wilson fermion as a function of quark chemical 
potential /j, and temperature T. We find evidence of a superfluid phase at intermediate /i and low 
T where the quark number density and diquark condensate are both very well described by a Fermi 
sphere of nearly-free quarks disrupted by a BCS condensate. Our results suggest that the quark con- 
tribution to the energy density is negative (and balanced by a positive gauge contribution), although 
this result is highly sensitive to details of the energy renormalisation. We also find evidence that 
the chiral condensate in this region vanishes in the massless limit. This region gives way to a region 
of deconfined quark matter at higher T and fi, with the deconfinement temperature, determined 
from the renormalised Polyakov loop, decreasing only very slowly with increasing chemical potential. 
The quark number susceptibility \q does not exhibit any qualitative change at the deconfinement 
transition. We argue that this is because \q ls n °t an appropriate measure of deconfinement for 
2-color QCD at high density. 



I. INTRODUCTION 

Despite over a decade of intensive efforts to unveil the 
phase structure of strongly interacting matter at high 
density (beyond a few times the nuclear saturation den- 
sity) and low temperature, even the question of which 
phases exist remains unanswered. A quantitative knowl- 
edge of this region would allow us to answer many ques- 
tions regarding the structure and properties of compact 
stars, including the question of whether deconfined quark 
matter can exist inside such stars. The reason for the lack 
of definite progress on this issue is that standard weak- 
coupling methods arc inapplicable except at asymptoti- 
cally high densities, while the various model approaches 
that have been employed have not been sufficiently con- 
strained by input from experiment or first-principles the- 
oretical calculations to yield reliable information in the 
region of interest. Thus, while a wealth of information 
exists regarding possible phases and their properties in 
various models, no reliable, quantitative results are avail- 
able as yet. For a recent review of high-density QCD, see 
Ref. [j. 

Many of the outstanding questions could in principle 
be answered by lattice QCD simulations, but these have 
been hindered by the notorious sign problem. While no 
method has as yet been shown to solve the sign problem 
for QCD, lattice simulations may still constrain model 
calculations by providing first-principles, nonperturba- 
tive results for QCD-like theories without a sign problem. 
This is the main aim of the present study. 

Among these theories, QCD with gauge group SU(2) 
(two-color QCD or QC2D) is of particular interest in 
that it shares most of the salient features of real QCD 
(eg, confinement, dynamical chiral symmetry breaking 
and long-range interactions). It differs from QCD in 
that the baryons of the theory are bosons, and the light- 
est baryon is a pseudo-Goldstone boson, degenerate with 
the pion (note though, that SU(2) models with adjoint 
matter Q and G2 with fundamental matter ||, both of 



which arc free from a sign problem, are expected to con- 
tain fcrmionic baryons in the physical spectrum). There- 
fore, instead of a normal nuclear matter phase this the- 
ory has a superfluid state characterised by condensation 
of these baryons, which at this point become true Gold- 
stone bosons. This has been observed in a number of 
lattice simulations; in particular, the excitation spec- 
trum including the Goldstone bosons has been studied in 
Refs. [J, [5| . A transition to a state of deconfined quark 
matter is expected at high chemical potential [i (see how- 
ever [H ) , and evidence of this was found in 0, H| . The 
precise nature of this transition remained unclear, how- 
ever, and in this paper we will attempt to answer some 
of the outstanding questions about this. 

An intriguing possibilty is that in an intermediate 
regime, strongly interacting matter may enter a chirally 
symmetric and confined phase, dubbed quarkyonic @. In 
3, it was suggested that the scaling of thermodynamic 
quantities with fx in the intermediate regime could be a 
sign of such a phase. It was not possible to draw any 
further conclusions, not least because the presence of a 
non-zero diquark source j 7^ 0, introduced to stabilise the 
simulations, distorted the /i-dcpendence of the relevant 
quantities. This will be remedied in the present paper. 

This paper is organised as follows. In Section [IT] we 
present results from simulations at zero chemical poten- 
tial. These results allow us to map out lines of constant 
physics, including the line of zero quark mass, which will 
in the future allow us to perform controlled extrapola- 
tions to the continuum and chiral limits, and also by 
varying N T at fixed cutoff to estimate the critical temper- 
ature Td for deconfinement at n = 0. In addition, these 
results form a large part of the input into the renormali- 
sation of energy densities, which is described and carried 
out in Section IIIII Section IIVI contains the bulk of our 
results for the (/u, T) phase diagram. In Section HV Al we 
present results for the order parameters for superfluid- 
ity and deconfinement, giving us an outline of the (//, T) 
phase diagram. Section IIVBI contains results for the 
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TABLE I: Simulation parameters, pi and rho masses and lat- 
tice spacing at li = j = 0. 



thermodynamic quantities, baryon density and (renor- 
malised) energy density, while Section IIV CI contains re- 
sults for the quark number susceptibility (preliminary re- 
sults from this work were presented in Ref. [HI), and 
in Section llVDI we investigate chiral symmetry breaking 
and restoration. Finally, in Section [V] we summarise our 
results and their implications. 
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TABLE II: Critical hopping parameter n c given by m 7r (K c 
0, for different values of /3. 



(pscudoscalar meson) mass , ratio of pion to rho (vec- 
tor meson) mass m^/mp and lattice spacing a. The lat- 
tice spacing was determined by fitting the static quark 
potential to the Cornell form V(r) = C + a/r + or and 
taking the string tension to be \fo = 440McV. 

We can determine the value k c ((3) where the quark 
mass vanishes by performing a linear extrapolation of 
m^. in 1/k for each value of /3. The results of this are 
shown in Table [TTJ 

We have also investigated the thermal deconfinemcnt 
transition at /z = using the fixed-scale approach. We 
have generated configurations with N T = 4 — 10 at j3 = 
1.9, k = 0.168, corresponding to a temperature range of 
113-281 MeV. At each temperature we have computed 
the Polyakov loop (L), which is an order parameter for 
dcconfincment of static color charges in the pure gauge 
theory, and exhibits a rapid crossover in a theory with 
dynamical fcrmions. It is related to the free energy F q of 
a static quark by 



L = e 



-F,(T)/T 



(3) 



II. SIMULATION DETAILS AND VACUUM 
PHASE STRUCTURE 

We study QC2D with a conventional Wilson action for 
the gauge fields and two flavours of Wilson fermion. The 
fermion action is augmented by a gauge- and iso-singlet 
diquark source term which serves the dual purpose of 
lifting the low-lying eigenvalues of the Dirac operator and 
allowing a controlled study of diquark condensation. The 
quark action is 

Sq +Sj=J2 ^ + [V>2 r (C75)T 2 Vi - h.c], (1) 



The free energy F q is only defined up to an additive 
renormalisation constant AF, which depends on the 
bare couplings /3, k. Different prescriptions for determin- 
ing this constant correspond to different renormalisation 
schemes. We have imposed the condition that the renor- 
malised Polyakov loop on our N T = 4 lattice (T = 263 
MeV) is equal to 1, or in other words, the free energy 
is zero at this temperature. We can then compute the 
renormalised Polyakov loop Lr(T) at any other temper- 
ature T from the bare Polyakov loop Lq via 



Lr(T) 



-F r {T)/T 



L (T)e 



-AF/T 



-(F a (T)+AF)/T 



Z^L (T= 1/aNr). 



(4) 



with 

M xy = 8 xy - kY^\0- ~ lu)e^°U^(x)S y ^ +i> 

+ (l + 70e"^"°t^(y)*»,*-4 (2) 

Further details about the action and the Hybrid Monte 
Carlo algorithm used can be found in Q ■ 

We have performed an extensive exploration of the pa- 
rameter space in the vacuum (T = (/, = j = 0) in the 
range /3 = 1.7 — 2.1. The parameters used are shown in 
Table U together with the values obtained for the pion 



where Z L = exp(-aAF) = L (N T = 4)" 1 / 4 (this pro- 
cedure was first outlined in Ref. [Hj]). The results are 
shown in Fig. U as a function of aT = 1 /N T . The red 
curve in Fig. [T] is a cubic spline interpolation between the 
data points. Taking the derivative of this (denoted by the 
green curve), we find the maximum at Ta = 0.193. If we 
instead use an Akima spline to interpolate, the maximum 
of the derivative appears at Ta = 0.183. Taking the cubic 
spline as our best estimate and conservatively estimating 
the uncertainty to be twice the difference between the 
Akima and cubic spline estimates, our result for the de- 
confinement temperature is Td(fi — 0) is T^a = 0.193(20) 
or T d = 217(23) MeV. 
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FIG. 1: The renormalised Polyakov loop La as a function of 
temperature T, for 16 3 x 7V T lattices at ft = 1.9, k = 0.168, /j, = 
j — 0. The red band is a cubic spline interpolation between 
the data points, and the green curve shows the derivative of 
the interpolation curve, divided by a factor of 10. 



III. RENORMALISATION OF ENERGY 
DENSITIES 

To determine the energy density, it is convenient to 
introduce different lattice spacings a s , a T in the space and 
time directions, with an anisotropy parameter £ = a s /a T . 
The energy density e is then given by [TU, sec. 5.4.1] 



e(T) 



1 dZ 



VdT- 



i 



dS 



NfN T a 3 s a T \ d£ 



(5) 



where we have used V = (N s a s ) 3 , T 1 = N T a T , and 



_d_ 



as_d_ 



(6) 



The partial derivatives must be taken with all other phys- 
ical parameters kept fixed. In our case, this means that 
the physical quark mass, and therefore the ratio /m p , 
is kept fixed. 

The anisotropic action S = Sg + Sq + Sj describing 
Nf = 2 Wilson quark flavors is given by 



Sn — — — 

G N r 



— ^ Re Tr U l3 (x) + 7 3 Rc Tr u «> (x) 



7s 



x,i<j 



(7) 



+ Kj2r(x)(D^) a (x) (8) 
Sj =KjYM> 2tr {*)Cl^ l {x) - ^(x)C l5 T 2 j 2tr (x)} , 

X 

(9) 



with 

(D^) a (x) = ( 7i - l)Ui(xW*(x + t) 

-( 7i + i)u}(x-i)i> a (x-i), (io) 

(D ^) a {x) = ( 7o - l)U (x)e^ a (x + 6) 

-(j + l)tf (x-0)e-^ a (x-6). (11) 

We also define 



7 S 



(12) 



The parameters j g and j q are the bare gluon and quark 
anisotropics, which in our formalism will be taken to be 
independent. 

Substituting these expressions into ^ (and dropping 
the \ as from all partial derivatives as it will be under- 
stood), we then readily derive 



£ g _ f ( N T a T 
T 4 ~ H N.a. 



dS G 



(ReTr^-) U" 1 ^ 



(RcTr[/ l0 ) 7 3 



dJ3_ 



(13) 



This coincides with the first part of Eq. (17) of Ref. [13j j . 
The terms in angled brackets are the average spatial and 
temporal plaquettes respectively, and the terms multi- 
plying them are what are usually known as the Karsch 
coefficients. In the weak coupling isotropic limit /3 — ¥ 
oo, 7 9 = 1 we have 



dlg_ dig 1 d£_ = _ 5/3 

d£ d£ ' d£ a da ' 

and we recover the expression used in 0, H[ : 



(14) 



9 

T 4 



N c 



[(Re Tr U i0 ) - (Rc Tr Uij)] . (15) 



The quark contribution to the energy density is given 



by 



N T a T 

N.a. 



dSQ 



(16) 



8k 



The terms in angled brackets are calculated using a 
stochastic estimator. Note a potentially useful identity 

WiipDoiP) +Kj2(i>DiiP) + = - Tr 1 = -4N c N f . 

(17) 

Note that we have taken explicit account of the minus 
sign associated with closed fermion loops in the defi- 
nition of the bilinear expectation values, ie. {tpTtp) = 



4 



— Tr(rM _1 ). It is therefore sufficient to evaluate the 
first and third terms on the LHS, enabling the second 
term, which enters into Eq. (fT5)l . to be estimated. In the 
isotropic limit 7 9 = £ = 1 this reduces to 



= N 4 

JI4 ±y T 



(ANfN c + m)K- x ^ - K^tfDoi/,) 



(18) 

In the weak coupling isotropic limit dn/dt; = 0, d^ q /d^ = 
1 and we recover 



(19) 



which coincides up to an overall sign with the expres- 
sion in 0, HI , where the fermion's Grassmann nature was 
ignored. 

Finally, the diquark contribution is given by 



ji4 



N. 



2N? fdj jdn\ 
= l^{dt + - K d-d {qq) (20) 

in the notation of Q- However, in the U(l)_B-symmetric 
limit j ' — > the second term inside the brackets vanishes, 
and since this limit is always found at j = for any 
anisotropy £, the first term also vanishes here. 
Similarly, the trace anomaly is given by 



Tan = e - 3p = 



dS 1 
'da. 



(21) 



With our anisotropic action the quark and gluon contri- 
butions are given by 




x 9B Q fry" 1 ' 



7g a S~ + P a 



da 



dp 
^a 



j3a 



da J 



(T, 



aujq 



1 

e 



dn dj q 
+ (V^oV) ( 1^ + m— 



(22) 



(23) 



However, in the isotropic limit, the bare anisotropics are 
always 1, and hence the derivatives dj 9tq /da vanish. We 
are then left with the standard expressions for the trace 
anomaly, 



(Tuu) g = -o^ J- (Re Tr U l3 + Re Tr U i0 ) 



(24) 
(25) 



Eqs. (|24l25p differ from the expressions used in 0, H[ by 
an overall factor f3 and an overall sign respectively; the 
resulting error is corrected in this paper. 



So, in order to evaluate the full energy density (ignor- 
ing j 0) from Eqs. (|13|16|) we need the following, which 
go into the definition of the "Karsch coefficients" : 



(26) 



dj3 dj g dn d^ q 
<9£ ' <9£ ' d£, ' dt, ' 

These are computed using the method presented in 
[T3I [l4l | . In addition to the bare anisotropies we define 
the physical anisotropics £ 9 = a s /a T as determined from 
gluonic observables such as the "sideways potential" [ll| , 
and £ g = a s /a T as determined from a meson dispersion 
relation. For a parameter set corresponding to a physi- 
cal system £ g and £ q should be equal, since otherwise a 
massless meson would not propagate at the correct speed 
of light; choosing the bare parameters to bring this about 
is a non-trivial tuning problem fl3l . [Lsj . In attempting 
to calculate the Karsch coefneents for the parameter set 
f3 = 1.9, k = 0.168, we do not attempt this tuning, but 
rather simulate unphysical ensembles with either j g or 
7 9 set to unity; the parameters are given in Table Mil In 
addition we use the isotropic ensembles given in Table [J 
For each ensemble we compute the ratio M = 



(rn^/mp) 1 , the lattice spacing a 



the gluon 



anisotropy £ 9 (from the sideways potential) and the quark 
anisotropy £ q (from the pion dispersion relation). The 
quark and gluon anisotropies are combined to form the 
average anisotropy £+ = | (£ g + £ 9 ) and the anisotropy 
mismatch £_ = £ g — £ q . Each of these quantities is fitted 
to a linear function in the bare parameters, 



a — ao 
a 

M - M 



Mo 



= aiA7 g 
= a 2 Aj g 

= a 3 Aj g 
= a 4 A7 9 



hA lq 
hAj q 



Cl A/3 
c 2 A/3- 

c 3 A/3- 
C4A/3. 



o^Ak . 

d^An . 
diAn . 



(27) 
(28) 

(29) 
(30) 



o 4 A 7? 

where ao and Mq are the values of a and M at the ref- 
erence point ft = 1.9, k = 0.168, 7 g = "/ q = 1, and Ax is 
the deviation of the bare parameter x from its value at 
the same reference point. Inverting the 4x4 matrix of 
coefficients (a,, bi,Ci,di) gives us the "generalised Karsch 
coefficients", which are the derivatives of the bare pa- 
rameters with respect to the "physical" parameters [31j 
£+, a, M. The first column gives us the Karsch coef- 
ficients (j2"6"|) . while the second column gives us the beta 
functions df3/da,dK,/da. 

Since we do not need to renormalise the pressure, 
knowledge of the beta-functions is not required here. 
However, we can use information about them to per- 
form consistency checks. In the isotropic limit, two of 
the Karsch coefficients can be expressed in terms of beta- 
functions, since 



dp_ 



£=1 



dfi_ 
da 



Oh 
~dj 



= -ap. (31) 
da 



We can also independently estimate the beta-functions 
from the isotropic results in Sec. [TH by taking derivatives 
wrt a along lines of constant physics. 
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Ps /3t Ks Kt 7s 7q 


£ g £ 9 m^jmp a s (fm) 


1.90 1.90 0.1680 0.1680 1.0 1.0 


0.968±l 1.07±1 0.807±1 0.178±6 


2.37 1.52 0.168 0.168 0.8 1.0 
1.27 2.83 0.168 0.168 1.5 1.0 
1.90 1.90 0.180 0.157 1.0 0.87 
1.90 1.90 0.147 0.192 1.0 1.3 


+ 2 +14 +4 +4 

0.720-2 0.853-io 0.805-s 0.177-3 
1.321±5 1.32±l 0.648-12 0.125-6 
0.747±4 0.78±3 0.746±?3 0.107±1 
1.146±4 1.53±1 0.946±! 0.229-13 


1.80 1.80 0.1740 0.1740 1.0 1.0 
1.90 1.90 0.1685 0.1685 1.0 1.0 
2.00 2.00 0.1620 0.1620 1.0 1.0 
2.00 2.00 0.1630 0.1630 1.0 1.0 


+3 +1 +6 +6 

0.989-3 1.03-3 0.777-8 0.177-8 
0.945-e 0.98-3 0.760±i8 0.153-is 
0.921±5 0.99±4 0.829±l 0.166±3 
0.881±5 1.04±I 0.773±i! 0.148±? 



TABLE III: Anisotropic lattice parameters and anisotropy results. The uncertainties are purely statistical. 




FIG. 2: The pion dispersion relation from the anisotropic 
12 3 x 24 lattices in Table ITTT1 



Results for the observables on our anisotropic lattices 
as well as the isotropic lattices used in this study, are 
given in Table [TTTJ Figs [2] and [3] illustrate the determi- 
nation of the quark and gluon anisotropics respecti vely . 
The gluon anisotropy in Fig. [3] was computed using [171 ] 

V xt (R 2 )-V xt (Rx) 



V xy (R 2 )-V xy (R 1 ) 



(32) 



where V x t (x) , V xy (x) are the potentials obtained from 
Wilson loops in the (x,t) and (x 7 y) plane respectively, 

W„(x,y) ~ Z xy e- yV ^ x \W st (x,t) ~ Z xte - tv ^ , 

(33) 

which is valid for large x and t,y. The fermion anisotropy 
is determined from the pion dispersion relation, 



2 zt2 2 2 



1 2 



(34) 



Hence, a straight-line fit of a 2 E 2 vs a 2 p 2 , as shown in 
Fig. [21 will give the anisotropy £ q . 

The results of the fits to (J27J) ([301) are shown in Ta- 
ble IIV1 We see that the \ 2 per degree of freedom is very 
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p\=P t = 1.9, 1^ = ^ = 0.168 

1.9,1c =0.180, K, = 0.157 
2.37,^ = 1.52.K =K=0.168 




FIG. 3: The gauge anisotropy from the anisotropic 12 x 24 
lattices in Table ITTT1 computed according to Eq. (I32[) , 
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TABLE IV: Results for the fits to Eqs. ((27 
the x 2 P er degree of freedom for each fit. 



M- X /N d f is 



high, especially for the average anisotropy and the mass 
ratio fits. This indicates that our linear approximation 
breaks down in this region, something which in the case 
of the anisotropy may be seen directly from the numbers 
in Table IIII1 where a nonlinear response of the physical 
anisotropics (and, indeed the lattice spacing) to the bare 
anisotropics is evident. To account for this, we would 
need to either include nonlinear terms in our Ansatz or 
employ smaller anisotropics (which would again require 
much higher statistics to determine the coefficients with 
sufficient precision). That is beyond the scope of this 
study. 

The generalised Karsch coefficients are presented in 
Table [V] We see that although the anisotropy deriva- 
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TABLE V: Results for the generalised Karsch coefficients 
dci/dxi. The numbers in the first column are the actual 
Karsch coefficients, while the second column gives the beta 
functions. 
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+ 17 
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TABLE VI: Results for the beta functions add /da and mass 
derivatives Mdci/dM , computed from fits to the isotropic 
data sets. 
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FIG. 4: The number of conjugate gradient iterations N cg 
per inversion, step size dt (solid lines) and acceptance rates 
(dashed lines) for our simulations on N T = 24 lattices. 



tives are reasonably well determined, other quantities, 
including the beta functions, have quite large uncertain- 
ties. The same has been found previously in real QCD 
with anisotropic lattices [l4| • It is likely that the extrac- 
tion of the lattice spacing from the static quark potential 
is the main limiting factor here, and that a high-precision 
lattice spacing determination from for example the Wil- 
son flow 18[ would help in this respect. 

A surprising result is the small value for the coefficient 
djg/d^, which comes out between 0.1 and 0.2, in contrast 
to c>7 g /<9£, which has a value close to 1 as expected. It is 
possible that this is related to the breakdown of the lin- 
ear approximation, and that including non-linear terms 
might bring this coefficient closer to 1. As we shall see in 
Sec. IIVB| this has a significant impact on the resulting 
energy density. 

The coefficients ad^ gA jda should be zero in the 
isotropic limit. While consistent or nearly consistent with 
zero within errors, the central values in Table fVl arc fairly 
large. If we could constrain these to be exactly zero, our 
overall uncertainties might be reduced. We also see that 
Eq. (|3"Tj) is satisfied within the admittedly large uncer- 
tainties. Again, it might improve the accuracy of our 
determination if these equations could be constrained to 
hold exactly. 

We may also use m T /m p instead of M = (m^/mp) 2 as 
our mass observables in the fits. We find that repeating 
the analysis above with this choice does not change the 
results for the Karsch coefficients and beta functions by 
much. 

We have also computed the beta functions separately 
from a 2-dimensional fit to the isotropic ensembles only. 
The results are shown in Table I VII As we can see, the 
two approaches give consistent results, suggesting that 
the systematic uncertainties of the method are under 
reasonable control. The numbers are also roughly consis- 



tent with (but somewhat larger than) the crude estimates 
used in Ref. [H, where a simple backward derivative ap- 
proximation was used. 



IV. RESULTS AT fi / 

We now focus on the (/? = 1.9, k = 0.1680) parameter 
set, and explore the interior of the (T, /i) plane for these 
bare couplings. Results for j = 0.04 on the 12 3 x 24 
lattices were already presented in Now, with the 
addition of data for j = 0.02 and, for some selected \x- 
values, j = 0.03, we can extrapolate all our results to the 
j = limit, assuming a linear behaviour with j. With 
some exceptions we do not see significant deviations from 
a linear behaviour in the cases where we have three j- 
values; the exceptions and their possible implications will 
be discussed below. We have also attempted to fit the 
data with a pure power law Ansatz, which describes the 
data well only in the region [i < \i w 0.3a -1 , and a 
power law plus constant form y = Aj a + B which does 
not yield any stable fits. 

We have also explored higher temperatures with data 
at N T — 16, 12, 8, and studied finite volume effects with 
the addition of a 16 3 spatial volume. The temperatures 
are T = 47, 70, 94 and 141 McV for N T = 24, 16, 12 
and 8 respectively. Details of our data sets are given in 
Tables IVIII Ixl Figure [4] shows the computational effort 
for the N T = 24 lattices in terms of the number of conju- 
gate gradient iterations per inversion and the molecular 
dynamics stepsize. It is evident from this figure that 
simulations in the dense region at the lowest j-value are 
1-2 orders of magnitude more costly than those of the 
vacuum. 
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0.95 


257 








1.00 


250 




600 




1.10 


252 




504 





TABLE VII: Number of trajectories for fj, 0, f3 = 1.9, n = 
0.168, N T = 24 (T = 47 MeV). The ja = 0.02, 0.03 configura- 
tions all have N s = 12. All trajectories have average length 
0.5. 
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0.900 
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TABLE VIII: Chemical potential values and number of tra- 
jectories for the 12 3 x 16 lattices (T = 70 MeV). The diquark 
source is ja = 0.04 in all cases. All trajectories have average 
length 0.5. 



afj, 

A(0.04) 



aft 

A(0.04) 
A(0.02) 



0.200 0.250 0.275 0.300 0.325 0.350 0.360 0.375 0.390 
1000 2500 2520 2520 2800 4900 2100 4900 1300 



0.400 0.425 0.450 0.500 0.600 0.700 0.800 0.900 
1080 1050 1050 1164 1128 600 510 540 
500 512 310 300 250 255 



TABLE IX: Chemical potential values and number of trajec- 
tories for the 16 3 x 12 lattices (T = 94 MeV). iV(0.04) and 
iV(0.02) are the number of trajectories for ja = 0.04 and 0.02 
respectively. All trajectories have average length 0.5. 
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TABLE X: As Table [IX] for the 16 3 x 8 lattices (T = 141 
MeV). 
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FIG. 5: The diquark condensate (gg)/// 2 extrapolated to j 
for N T = 24, 12, 8 (T = 47, 94, 141 MeV). 



A. Order parameters and phase structure 

Figure [5] shows the diquark condensate, 

(gg) = (ip 2tr C 75 T 2 ^ - ^C l5 T 2 j 2tr ) - h.c. , (35) 

as a function of chemical potential, for the N T = 24, 12 
and 8 lattices. In the case of a weakly-coupled BCS con- 
densate at the Fermi surface, the diquark condensate, 
which is the number density of Cooper pairs, should 
be proportional to the area of the Fermi surface, ie 
(gg) ~ [i 2 . This is to be contrasted with chiral pertur- 
bation theory (%PT) [l9| . which for fx \i Q at leading 
order predicts (gg) to be //-independent. 

For the lowest temperature T = 47 MeV {N T = 24) 
we see an almost perfect proportionality in the region 
0.35 < /ia < 0.6. The lower limit of this region roughly 
coincides with the onset chemical potential /i D m 
0.33a _1 , below which both the quark number density and 
diquark condensate is expected to be zero. The reason 
we see a gradual rise from [ia ~ 0.25 is our use of a linear 
Ansatz for the j — > extrapolation, which is not valid 
in this regime. We have confirmed that a pure power 
law, (gg) = Aj a , works satisfactorily here. The power 
a is found to decrease monotonically with /z, from 0.86 
at fia = 0.25 to 0.56 at fia = 0.35. At fia = 0.30, the 
only point in this region where we have 3 j-values at 
our disposal, we find a = 0.709(6). At fia = 0.5, the 
X 2 for the linear fit to 3 j-values is 3.3, while a power- 
law fit gives x 2 = 2.0, suggesting that a power-law plus 
a constant might be preferable. However, this does not 
lead to any significant change in the extrapolated value, 
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FIG. 6: The renormalised Polyakov loop as a function of 
chemical potential, for all temperatures. The open symbols 
are for ja — 0.04; the filled symbols are extrapolated to j = 0. 
The inset shows the unrenormalised Polyakov loop. 
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FIG. 7: A tentative phase diagram, including the location of 
the deconfinement transition in the (/i, T) plane, determined 
from the renormalised Polyakov loop, and the transition to 
the diquark condensed (qq) ^ phase. Also shown is the 
deconfinement point from Ref. 0- 



and hence we choose to quote our results from the linear 
fit also here. 

For \xa > 0.6, (qq) / \i 2 rises again before possibly reach- 
ing a new plateau at /js« 1.0. This is possible evidence 
of a transition to a new state of matter at high density, 
but at these large densities the impact of lattice artifacts 
cannot be excluded. 

At T = 70 MeV (N T = 16) we are not in a position to 
perform a j — > extrapolation, but from the ja = 0.04 
data we see only a mild suppression in (qq) , and only for 
Lia > 0.8. Since the results are almost indistinguishable 
from those at T = 47 MeV we do not show them here. 

At T = 94 MeV (N T = 12) we see that (qq) is signifi- 
cantly smaller for all values of [i and drops dramatically 
above \ia > 0.7. This gives us the first indications of the 
transition between the diquark-condensed and the nor- 
mal phase. At T = 141 MeV (N T = 8) we find that 
the diquark condensate is zero at all [i, confirming that 
the system is in the normal phase at this temperature. 
A systematic investigation including more temperatures 
and an extrapolation to j = at all temperatures will 
be required to establish the exact location and nature of 
this transition. 

Finally, comparing the numbers from the 12 3 x 24 and 
16 3 x 24 lattices, no evidence of any significant finite vol- 
ume effects are found, except at [ia = 0.9 where the 
condensate on the smaller volume is slightly suppressed. 

Figure [6] shows the order parameter for deconfinement, 
the Polyakov loop (L), for our four different tempera- 
tures. It has been renormalised using ([4]), using the [i- 
independent rcnormalisation constant Zl already com- 
puted in SecHIl We see that for each temperature T, (L) 
increases rapidly from zero above a chemical potential 
fid(T) which we may identify with the chemical potential 
for deconfinement. However, since L is a convex function 
of fj, at all T, it is not possible to use the variation of 
L with fi to define /id(T). In the absence of a more rig- 



orous criterion, we have taken the point where L crosses 
the value it takes at = 0), Ld = 0.6, to define Hd(T). 
The results are shown in Fig. [3 with error bars denoting 
the range Ld =0.5-0.7. To more accurately locate the de- 
confinement line, wc will need to perform a temperature 
scan for fixed /a- values, as was done for fi = 0. 

For our lowest temperature (N T = 24), the renor- 
malised Polyakov loop is too noisy for any quantitative 
conclusions to be drawn. This is because the signal 
(which is consistent with for [ia < 0.75) as well as 
the statistical noise are multiplied by the large renormal- 
isation factor Z£ 4 = 2084. However, the unrenormalised 
Polyakov loop Lq, shown in the inset of Fig. [6l exhibits 
the same qualitative behaviour as for the higher temper- 
atures. We also find that there are no significant volume 
effects, while the diquark source tends to suppress the 
Polyakov loop slightly. At /ia w 0.75 we see that the 
curves for the renormalised Polyakov loop at the differ- 
ent temperatures cross, so that at higher fj,, L is smaller 
for higher temperatures. This, however, depends on the 
rcnormalisation scheme: if we had instead imposed the 
condition that Lr = 0.5 at N T — 4,/i = 0, the curves 
would not cross. 

The estimates of critical chemical potentials for both 
deconfinement and superfluidity can be translated into a 
tentative phase diagram, shown in Fig. [7] It is worth reit- 
erating that the points on the phase boundaries are rough 
estimates only, since we do not have a precise criterion 
for the transition. In Section IIV CI we will present re- 
sults for a different measure of deconfinement, the quark 
number susceptibility. Wc also show the estimate from 
the coarser lattice in Ref. 0. Clearly, a combination of 
temperature effects and lattice artefacts is responsible for 
the discrepancy between the /Zd-values quoted in @, @] ■ 

In Fig. [7] we also show our estimate of the transition 
between the supcrfluid and the normal phase. Again, 



9 



2.5 



3 « 



1.5 



0.5 



5 ca 

u a" 2 



O N=24 
♦ N =12 



A N. 



2_d 



in 



± 



3 



0.2 



0.4 



0.6 0.8 



FIG. 8: The quark number density at j — 0, divided by the 
density for a noninteracting gas of lattice quarks (top) and 
continuum quarks (bottom). 



since we do not yet have j ~ > extrapolated data at 
all temperatures, and because our temperature grid is 
fairly coarse, these transition points are also only rough 
estimates. 

In summary, from the order parameters we find sig- 
natures of three different regions (or phases): a nor- 
mal (hadronic) phase with (qq) = 0, (L) rj 0; a BCS 
(quarkyonic) region with (qq) ~ fi 2 at low T and inter- 
mediate to large /i; and a deconfined, normal phase with 
(qq) = 0, (L) 7^ at large T and/or /i. We cannot exclude 
a deconfined supcrfluid phase with (L) > 0, (qq) ^ at 
large ^ and intermediate T. 

After extrapolating our results to zero diquark source, 
we see no evide nce of a B EC region described by xPT, 
with (qq) ~ ~ Jif/J? in contrast with earlier 
work with staggered lattice fermions Q ■ This may be be- 
cause we do not have a clear separation of scales between 
the Goldstone diquark scale and more massive states, 
and hence the region of tightly bound diquarks is very 
narrow. A more pessimistic scenario is that the BEC 
region is masked by the poor chiral properties of Wil- 
son fermions. Simulations with lighter quarks may help 
clarify this. 



B. Equation of state 

We now turn to the bulk thermodynamics of the sys- 
tem: the quark number density n q , the pressure p and 
the energy density e. Figure [5] shows the quark num- 
ber density n q for N T = 24, 12 and 8, extrapolated to 
zero diquark source. The linear extrapolation is marginal 
(x 2 = 3.7) at N T = 24, /xa = 0.5; however, a quadratic 
form n q = uq + cj 2 does not work much better, and a 
constant + power law does not yield a stable fit. With- 
out higher statistics and additional j-values we are not 
in a position to determine the correct j-dependcnce here, 
and hence our choice is to use the linear extrapolation 
throughout. 

In the top panel of Fig. (8) we have normalised the den- 
sity by the density nigg for noninteracting fermions on the 
same lattice volumes (12 3 x 24, 16 3 x 12, 16 3 x 8). This 
is the same normalisation as was used in @, H| . In the 
bottom panel, we have instead divided by the continuum, 
infinite- volume expression for noninteracting fermions at 
the same temperature and chemical potential. The dif- 
ference between the two gives an indication of the lat- 
tice artefacts. We see that the density rises from zero at 
/i rj fi = 0.32a -1 , and for the two lower temperatures is 
roughly constant and approximately equal to the nonin- 
teracting fcrmion density in the region 0.4 < fia < 0.7. 
The peak at ijlcl ~ 0.4 in the N T = 24 data in the up- 
per panel is an artefact of the normalisation with usb 
for a finite lattice volume. This is evident from Fig. [9l 
where in the upper panel we show the quark number 
density at ja = 0.04 for different lattice volumes, nor- 
malised by nsB for the corresponding lattice volumes, 
while in the lower panel we have used the same normal- 
isation for all lattices. We have chosen to use n$B for 
a 16 3 x 24 lattice for this normalisation; note that this 
choice is purely a matter of convenience, the purpose be- 
ing to easily compare the raw numbers for n q from differ- 
ent lattices. We see that there is no difference between 
our raw numbers for n q on the 12 3 x 24 and 16 3 x 24 
lattices at j = 0.04; however usb for the 12 3 lattice has 
a dip around fia ~ 0.4, which is absent for the 16 3 lattice 
(see below). We therefore conclude that our previous 
interpretation Q of the peak in n q JnsB in this region 
as evidence of a BEC condensate described by xPT was 
probably erroneous. 

The density for N T = 8 does not show any plateau as 
a function of /x; instead, n q /nsB shows a roughly linear 
increase in the region 0.4 < /ia < 0.7. This is suggestive 
of the system being in a different phase at this tempera- 
ture. We also note that n q /nsB for N T = 12 rises above 
the corresponding N T = 24 data for //a > 0.7, where, 
according to the results of Sec. IIV Al the hotter system 
is entering the deconfined, normal phase. 

Our results for N T = 16 are indistinguishable from the 
N T = 24 results except for fia > 0.8, where they also 
start increasing above the N T = 24 values. The rise in 
n q/ n SB for A* a ^ O-^ rnay be a signal of a new phase, 
although in this region the influence of lattice artefacts 
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FIG. 9: The quark number density at ja = 0.04 for different 
lattice volumes, divided by the density for a noninteracting 
gas of lattice quarks on the same volume (top) and on a fixed 
volume of 16 3 x 24 (bottom). 



cannot yet be ruled out. 

These results lend further support to our previous con- 
jecture that in the intermediate-density region the system 
is in a "quarkyonic" phase: a confined phase (all excita- 
tions are colourless) that can be described by quark de- 
grees of freedom. We reiterate that because of the large 
explicit breaking of chiral symmetry in our simulations, 
we cannot say anything at this point about chiral symme- 
try restoration, another characteristic of the quarkyonic 
phase conjectured in Ref. We will come back to this 
issue in Sec. IIVDI 

Further insight on the influence of both UV and 
IR artefacts can be gleaned by considering the ratio 
n ss /^Sff = 0)' calculated for two different volumes 
using the formula given in Q , and shown in Fig. 1101 It 
should be noted that the correction is numerically big 
across large portions of the /i-axis. The oscillatory be- 
haviour seen for fia < 0.8 is an IR artefact known to 
arise from the non-sphericity of the Fermi surface result- 
ing from the discretisation of momentum space [2fjj |. In 
particular, the dip at fj,a ~ 0.4 on 12 3 x 24 coincides with 
the peak in n q /nsB seen in the upper panel of Fig. 
suggesting that it is an artefact of the attempted lattice 
correction; indeed, the peak is not present on 16 3 x 24 
since the dip in the correction factor has moved to smaller 
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FIG. 10: Ratio n^/n^™* evaluated for free massless quarks 
on both 12 3 x 24 and 16 3 x 24 lattices. The inset shows the 
same ratio for the 12 3 x 24 lattice, for four different values of 
the diquark source j. 



H . By contrast, the correction factor coincides on the two 
volumes for fia > 0(1), suggesting that the considerable 
departure from unity at large fj, is due to UV effects. As 
we can see in the inset of Fig. [101 the diquark source has 
a negligible effect on the noninteracting quark density, 
and hence the strong j-dcpendcncc we sec in our results 
arises from the interactions. 

Next we discuss pressure, which as the negative of the 
free energy density, may be calculated via the integral 
of any thermodynamic observable along an appropriate 
contour. It is particularly convenient to integrate along 
the ^-axis via p = f 1 n q dn, since the cutoff does not 
change. Here hq is chosen so that p(p.o) = to good 
approximation; in the limit T — > /io should coincide 
with the onset fi . 

In our analysis the integral is readily approximated by 
a trapezoidal rule; as always, we present data normalised 
by the free field value psb, a procedure not uniquely de- 
fined away from the continuum limit. We have examined 
three schemes: 



- 

\PsbJo 

(-) 

\PSB J j 

(—) 

\PsbJu 
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J Ho 
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(psB^r 1 



H cont 
n SB 
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(^')nq(fJ,')dfj,' . 



(36) 
(37) 

/ 

) 

(38) 



where 



..cont NfNc f 4 „ q ,,2rp2 , ^_rp4 \ (or\\ 

PsB = 12^ ^ T + j (39) 

is the continuum pressure for a free gas of quarks, and 
p l gg the corresponding value obtained by summing over 
free quark modes on the finite lattice. Versions (|37p and 
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FIG. 11: p/psb vs. £ta for ja = 0.04 and various temperatures. Also shown are values extrapolated to j = for N T = 24. Top 
left: (p/Psb)o,i,ii for j -)■ 0, N T = 24. Top right: (p/psb)o- Bottom left: (p/psb)i- Bottom right: (p/psb)ii- 



(|3"8")) were both studied in [tJ, whereas only (p/psb)ii was 
used in [|,[2l| . 

Fig. [11] for data taken with ja — 0.04, as well as the 
j —> extrapolated data for N T = 24. In scheme II 
lattice data are "corrected" for artefacts before integrat- 
ing. The results clearly inherit the bump at fia ~ 0.45 
also manifest in Fig. [9j which we now believe to be an 
IR artefact. However, this bump is absent (or strongly 
suppressed) in the j — > limit, mirroring the absence 
of a significant bump in the upper panel of Fig. [8] By 
contrast the scheme data have the ratio p/psB sub- 
stantially exceeding unity in the large- \x regime above 
deconfinement, which probably reflects the fact that UV 
artefacts are not being fully corrected here. For this rea- 
son we now prefer scheme I, where for the coldest lattice 
p/psb has a plateau with value ~ 1.4 in the suspected 
quarkyonic region, only rising to ~ 2 at large fi. Again, 
therefore, we conclude that for low T there is a range of 
fi where thermodynamic quantities scale approximately 
the same as free quarks; that the evidence for a peak 
above onset matching the expectations of xPT has sub- 
stantially diminished; and that p/psb rises above unity 
in the deconfined regime. By T = 141 MeV (N T = 8), 
however, the ratio rises monotonically and the distinction 
between these different regimes is largely washed out. It 



is clear, however, that the full story will only emerge once 
the continuum and thermodynamic limits arc both taken 
with care. 

The quark and gluon contributions to the energy den- 
sity, for ja = 0.04, are shown in the upper panel of 
Fig. [I2j We see that the quark energy density is al- 
most independent of temperature for all temperatures, 
while the gluon energy density shows a clearly different 
behaviour only for the highest temperature. We find that 
the gluon energy density is independent of the diquark 
source within errors, so these results are representative 
for the j — > extrapolated data. The quark contribution 
is sensitive to the diquark source in the low-/i region, as 
can be seen from the j : — > extrapolated data also shown 
in Fig. [H 

Comparing these results with the unrenormalised (and 
unextrapolated) results in Figs 1 and 3 of Ref. Q , we see 
a dramatic difference. Clearly, the proper renormalisa- 
tion is crucial to any reliable determination of the energy 
density, and in particular it is clear that the terms pro- 
portional to d(3/d£, and dn/dt; in (fT3"| and (fll)j) respec- 
tively cannot be ignored. To illustrate this more clearly, 
we show in Fig. [T3]the quark contribution to the energy 
density on the 12 3 x 24 lattice at ja = 0.04, computed us- 
ing different values for the Karsch coefficients. The open 
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various temperatures, for ja — 0.04 (open symbols) and ex- 
trapolated to j = (filled symbols). Bottom: Total energy 
density divided by /x 4 . 



circles correspond to the unrenormalised energy density 
which was presented in Ref. @ (note that the normal- 
isation is different). The other data sets correspond to 
different values of dj q /d^, with dn/d^ set to the value of 
—0.052 that was determined in Sec. IIII1 We have chosen 
to use the tree- level value of 1, the value 0.131 deter- 
mined in Sec. IIII1 and a value of 0.8, which is similar to 
the value found for dy g /d£, and at the margins of our 
95% confidence interval. We see that using the correct 
(non-zero) value for <9k/<9£ is most important at low /x, 
where this alone changes the sign of e q . At large /x, the 
djq/dt; term will dominate, as it does at tree level. 

It should be noted that the uncertainties in the Karsch 
coefficients are not included in the total uncertainties in 
the plots shown here. On the basis of Fig. [13] one may 
conclude that these uncertainties will have an effect of 
O(100%) in the energy density. 

Although e q appears to be negative at least for low /x, 
and possibly for all /x-values considered here, the total 
energy density e = e g + e qi shown in the bottom panel of 
Fig. I12[ remains positive or consistent with zero every- 
where in the j — > limit. Although on the face of it a 
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FIG. 13: Renormalised quark energy density divided by /x 4 
density at T = 47 MeV (N T = 24), ja = 0.04, for different 
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negative value for e q is surprising, it is notable that the 
renormalised quark energy density shown in Fig. [12] has 
a qualitative resemblance to the unrenormalised energy 
density (fl"o| measured for QC 2 D with Nf — 4 Wilson 
quark flavors in Fig. 5 of Ref. [21(. The parameters used 
in that study correspond to a much finer lattice, with 
a I yfa having a value approximately one-third that used 
here. It is conceivable, therefore, that the Karsch coeffi- 
cients for Nf = 4 fall far closer to their free-field values, 
and hence their neglect in [2l[ is much better justified, 
reinforcing the conclusion that £ 9 (/x) < 0. 

Finally, we consider the trace anomaly, computed ac- 
cording to Eqs ([24]) . ([25]) . which is shown in Fig. [T4l 
With the correct expression (|25|) . we now find the quark 
contribution to be negative for all /x, whereas in 0, H| 
it had erroneously been presented as positive. Since the 
beta-functions only enter into the expressions as overall 
constants, and our updated values are not dramatically 
different from those used in Ref. 0, the qualitative be- 
haviour of the N T — 24, ja = 0.04 data is the same as 
previously reported in Ref. Q , apart from the sign of the 
quark contribution. 

For small and intermediate /x, the gluon and quark con- 
tributions have opposite signs and similar magnitudes, 
leading to a nearly vanishing total trace anomaly in the 
region OZeq/xa < 0.7. The gluon contribution decreases 
for /x > 0.5 and becomes negative for /xa > 0.75, while 
the quark contribution has a plateau for 0.5 < /xa < 0.75 
and increases rapidly in magnitude thereafter. This leads 
to a negative total trace anomaly at large /x, which cor- 
responds to the positive and increasing pressure p = 
- r w )/ 3 observed in Fig. [TO 

We see no difference in the trace anomaly between the 
two lowest temperatures, T = 47 and 70 MeV. At T = 94 
MeV and 141 MeV (N T = 12 and 8) the gluon contri- 
bution becomes larger (or less negative) and the quark 
contribution becomes more negative at large /x. The net 
effect of this, however, is to leave the total trace anomaly 
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FIG. 14: Top: gluon (shaded symbols, dotted lines) and 
quark (open symbols) contributions to the trace anomaly, at 
ja — 0.04. The filled symbols denote the quark contribu- 
tions extrapolated to j = 0. Bottom: Total trace anomaly 
divided by fi 4 , for ja — 0.04 (open symbols, dashed lines) and 
extrapolated to j = (filled symbols). 



nearly unchanged. 

We find that the trace anomaly depends only weakly on 
the diquark source for nearly all T and /j,. The main effect 
is to increase the gluon contribution at large [i and T, and 
to decrease the magnitued of the quark contribution at 
low T, for large and small (i. It is quite striking that 
there appears to be little or no dependence on cither 
temperature or diquark source in the region \i a < /i < 
0.55/c. 

Once again, it is instructive to compare with the 
Nf = 4 study of Ref. [2l|. In that case (see figs. 6 and 
8 of [2l| ) , after taking into account the incorrect sign for 
the quark contribution, the gluon unrenormalised contri- 
bution to is negative for all (i < /i^, while the quark 
contribution is positive, which is the opposite of what we 



observe here. However, this still leaves open the possibil- 
ity of the two contributions nearly cancelling, giving rise 
to nearly-conformal matter in the quarkyonic region. 
C. Quark number susceptibility 

In a mathematical sense the Polyakov loop is a well- 
defined signal for deconfinement, at least in pure gauge 
theories; physically it reveals something about the be- 
haviour of static color sources, which are well approxi- 
mated by heavy quarks, in a baryonic medium. Recent 
studies of a non-relativistic formulation of QC2D [22j 
took the first step beyond the static approximation, and 
revealed a non-trivial T- and /x-depedence for s-wave 
states formed from heavy quarks. Another observable 
related to confinement is the quark number susceptibil- 
ity Xq = dn q /dfj.. This observable is usually thought of as 
encoding the fluctuations in the baryon (or quark) num- 
ber, and is of particular interest as a measure of confine- 
ment or deconfinement of light quark degrees of freedom 
[23M26I ] . If quarks are confined inside hadrons, the fluctu- 
ations of the quark number and hence the susceptibility 
will be suppressed, since increasing the quark number en- 
tails exciting a baryon, which requires a large amount of 
energy. If quarks are not confined, it is possible to excite 
a single quark, which requires much less energy, giving a 
larger quark number susceptibility. 

This link between Xq and deconfinement is clear in the 
case of QCD, where all baryons are heavy. In the case of 
QC2D the situation is less clear, since the lightest baryons 
are the pseudo-Goldstone diquarks, and large fluctua- 
tions are possible even in the confined phase. Nonethe- 
less, it is of great interest to study fluctuations in quark 
number at large density and low temperature. The 
only previous such study is Ref. [27|], where the Dyson- 
Schwinger equation in the rainbow approximation was 
employed. Hence QC2D offers an opportunity for a first 
systematic non-perturbative study of Xq m this regime. 

For an ideal gas of massless (continuum) quarks and 
gluons, at temperature T and chemical potential /1, we 
have: 



l SB 



N f N c 



cont - N f N c 



XSB 



3 

H 

3 



3tt 2 
,2 



(40) 
(41) 



Now consider the quark action (JTJ) rewritten in the form 
^M^ : where we have introduced the bispinors ^ = 
(V>i,C- 1 T 2 ^ r ) tr , * = (^i,-^ r CT 2 ); see Ref. ^ for 
details. From the definition of Xq w ^ have: 



Xq 



dn q 
dfjL 
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*— — * 
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O/J, 



-d 2 M ■ 



(42) 
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From this equation we can identify four different terms: 
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(43) 
(44) 
(45) 
(46) 



The second term of Eq. (|42p yields two terms, Ti and 
Ci, because there are two ways to contract the spinors. 

The calculation of the traces is done using unbiased 
estimators, introducing N n complex noise vectors r\ with 
the properties: (r? x ) = and (r/ x Vy) = $xy For example, 
the determination of the trace, used for T\ and T2, is 
based on the following relation: 



Tr 



M 



_xdM 



dii 



N, 



( dM 



■My/3j;z~/k r lz7k 



xai;y/3j 



(47) 

Because two independent source vectors are required to 
compute T2 , we refer to this term as "disconnected" ; the 
other three "connected" terms need only one source vec- 
tor. 

It turns out that the connected term gives an impor- 
tant contribution to \ q at low and high values of the 
chemical potential and therefore cannot be considered 
negligible. Moreover, it changes sign around /ia ss 0.66. 
On the other hand, the terms T\ and T2 are equal within 
errors but with opposite sign, i.e., their net contribution 
is consistent with zero everywhere. 

In Fig. [15] we plot the ratio Xg/XssS f° r f° ur different 
temperatures, versus the chemical potential. For an ideal 
gas of quarks and gluons this ratio would be a constant, 
see Eq. (|4ip . and we see that an approximate plateau 
is actually present for a/i < 0.55, at least for the three 
lowest temperature; after this value we can see a sharp 
increase of \ q . The value of the plateau is Xq/XSB ~ 1-6 
which is higher than the ideal value of 1.0. Moreover, 
it is evident from this plot that \ q is T-indcpendent at 
low temperature, since there is no significant deviation 
in the behaviour of the three curves. This is to be con- 
trasted with the Polyakov loop in Fig. [6l which shows 
deconfinement for three different values of 11 = /id(T), as 
the temperature is varied. Only for the highest tempera- 
ture do we see a different behavour signalling a different 
phase. 

It is also instructive to compare the numerical results 
with the equations corresponding to Eq. (|4*Tj) but taking 
in account the finite volume and the lattice discretisation. 
In Eq. (26) of Ref. 0] , the expression for the quark num- 
ber density n 1 ^ for free Wilson fermions on the lattice 
is presented, from which Xsb ^ s easily obtained. Fig. [TH] 
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FIG. 15: Ratio Xq/x'sb versus fj,, for ja = 0.04. The vertical 
dashed line marks the position of fj, a . 



plots the ratio Xq/x^SB f° r two values of the quark mass 
used in the determination of the subtracted bare 

quark mass m q = 1 /2k — l/2n c and the 'constituent' 
quark mass m c = m p /2. In this case we observe a dif- 
ferent behaviour for a/i < 0.45, followed by a flat region 
compatible with a ratio of one, ie. the system is behaving 
as free fermions, with again an increase for higher values 
of /i. Clearly, in Fig.[16]thc value of the mass used for the 
free fermions has a quantitative effect, in that the value of 
the plateau is shifted when the mass is increased, but this 
does not change the qualitative considerations. The re- 
sults using m q are almost identical those obtained setting 
m = 0. These plots again confirm the above scenario: we 
do not see any abrupt change for Xq as a function of T, 
whereas the Polyakov loop becomes different from zero 
at a T-dependent lid- Again, something different seems 
to emerge for the highest temperature, T = 141 MeV 
(N T = 8). 

The effect of the diquark source is illustrated in Fig.fTTl 
where we show X 9 /x5s( TO c) for the 12 3 x 24 lattice and 
ja = 0.04, 0.02 and 0. We find that the diquark source 
only has a significant effect for low jj,, where it increases 
the value of Xq slightly. 
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FIG. 16: The ratio between the measured quark number sus- 
ceptibility at ja = 0.04 and the ideal value for lattice free 
fermions for two values of the fermion mass: m q = 0.05 (top) 
and m c = 0.42 (bottom). The vertical dashed lines mark the 
position of fi . 



D. Chiral Symmetry Breaking 

An important issue which we have been hitherto un- 
able to address is the chiral properties of the ground state 
once fj, > [i Q . The original description of the quarkyonic 
phase in SU(iV c ) gauge theory was in terms of a chi- 
rally symmetric but confined medium, ie. one in which 
the chiral condensate ('tpip) — > as the bare quark mass 
m — > Later this picture was modified; chiral sym- 
metry breaking via a translationally non-invariant "chiral 
spiral" was postulated in (28|. 

It is clearly desirable to determine the fate of chi- 
ral symmetry breaking for the case of QC2D by a lat- 
tice calculation. Indeed, (ip"^} was examined in early 
studies such as @ using staggered lattice fermions, and 
reasonable quantitive agreement found over a decade of 
quark mass with the prediction of leading order xPT for 
T — > 0, namely that for fi < /i G the chiral condensate is 
/i-independent, and for fi > fi a 



('(pip) (X 



(48) 




Unfortunately, since the global symmetries of stag- 



FIG. 17: The ratio between the measured quark number sus- 
ceptibility at different diquark sources j and the ideal value 
for lattice free fermions with a fermion mass of m c = 0.42. 
The vertical dashed lines mark the position of fi . 



gered fermions do not coincide with those of continuum 
QC 2 D I , these results are not directly applicable. In any 
case, no attempt was mde to explore beyond the regime 
of applicabiliity of xPT. 

We should note that in QC2D, chiral symmetry is bro- 
ken by (ipip) , but remains unbroken by the diquark con- 
densate (qq) defined in (pJ5"|) . In real (3-color) QCD, differ- 
ent patterns of diquark condensation may lead to color 
superconducting phases which are either chirally sym- 
metric or chirally broken, and in the latter case the chi- 
ral symmetry breaking may be signalled by (ipip) 7= or 
by a different order parameter. In the following, we will 
ignore these complications and focus exclusively on the 
chiral condensate 

However, our use of Wilson fermions precludes any di- 
rect study in the current simulation, since this formu- 
lation violates chiral symmetry explicitly. Our strategy 
therefore is to calculate a chiral order parameter using 
a fermion formulation with manifest chiral and baryon 
number symmetries using the gauge backgrounds ensem- 
bles generated with Wilson quarks. The disparity be- 
tween valence and sea quarks violates unitarity; we miti- 
gate this uncontrolled approximation by tuning the mass 
of the valence quarks so that the pion mass coincides with 
that used in the simulation; once fi 7= the onset tran- 
sition of the valence quarks should then at least coincide 
with the true value. 

Rather than the obvious choice of staggered fermions 
for the valence quarks, we found it expedient to use 
the existing code for Nf = 2 Wilson fermions with 
the parameter r (which has the conventional value of 
unity in ([2])) set to zero. For j = this is equiv- 
alent to eight identical staggered fermions with mass 
m = (2k) _1 . For non-zero lattice spacing and [i 7= the 
action has a U(8)<g>U(8) global symmetry which is bro- 
ken by m 7= (explicitly) or (ipip) 7^ (spontaneously) 
to U(8)y (the subscript denotes vectorlike), which incor- 
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porates the U(l)s of baryon number. A diquark source 
j ^ breaks U(8)y to a SU(2)<g)SU(2) which preserves 
isospin but no longer includes \J(1)b- 

By studying effective mass plots as the valence Ky was 
varied we found that Ky = 8.0 gave the closest match to 
the value m^a = 0.66(2) found for /3 = 1.9, k = 0.168. 
Fig. [TH] then shows the resulting chiral condensate as a 
function of /i for the various lattices studied. Note that 
ja = 0.04 throughout, since this was found to yield a less 
noisy and more stable signal - hence these results are not 
reproducible using pure staggered fcrmions. Two things 
are apparent; first the shape of the curve is in qualitative 
agreement with the old staggered results of @ over the 
whole range of /j, studied, and thus consistent with ([48]) 
assuming an onset /x D a ~ 0.3. Secondly, the results are 
independent of temperature even up to T = 141MeV 
(N T =8). 

It appears that the chiral symmetry properties of the 
dense phase are well-described by x?T. In a sense the 
issuse of "chiral symmetry restoration" in QC2D is aca- 
demic, since the onset scale is set on the assumption 
that chiral symmetry is explicitly broken. Nonethe- 
less, we can characterise the dense phase by whether 
\im mv ^ ('tp' l P( m v)) vanishes or not. We determine this 
by using three different values Ky == 8, 16 and 40, and 
observing that with the field normalisations implicit in 



Kj(tp^)i _ m 2 (qq) 1 1 = 1 (ipi>)o = 0; 

iij(#)a rn 1 (qq) 2 [< 1 (#> ^ 0, m 2 < mi. 



(49) 

Here qq denotes the scalar quark bilinear with conven- 
tional normalisation and (ipip)o is the chiral condensate 
in the massless limit. Fig. [lUl shows this ratio plotted for 
both (8,16) and (8,40) valence mass pairs on the 16 3 x 24 
lattice as a function of /1, and clearly indicates symme- 
try restoration for fx > fj, . Very similar plots are found 
for the other temperatures explored. We therefore con- 
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elude that the gauge field backgrounds at high baryon 
density in QC2D are consistent with chiral symmetry be- 
ing unbroken by a scalar condensate, alt houg h the exotic 
translationally-non-invariant scenario of [28j is not ruled 
out. 



V. CONCLUSIONS AND OUTLOOK 

We have carried out the first extensive exploration 
of the phase diagram of two-color QCD (QC2D) in the 
(T, fi) plane using first-principles lattice simulations. Our 
main findings are summarised in the tentative phase dia- 
gram of Fig. [7] We find evidence of three distinct regions: 

1. A vacuum/hadronic phase, with (qq) = 0, (L) ~ 
0, (ifiip) ^ 0, n q ss 0, at low T and \x < no = itt-^/2; 

2. A quarkyonic phase at low T and intermediate to 
large fj,, which is confined ((L) ps 0) and charac- 
terised by a chiral condensate which vanishes in 
the chiral limit, Stefan-Boltzmann scaling of bulk 
thermodynamic quantities (including a nearly van- 
ishing trace anomaly) and BCS scaling of the di- 
quark condensate; 

3. A deconfined quark-gluon plasma phase at high T 
(and/or large /i). 

The main difference from our previous studies is that 
the BEC region has disappeared as a consequence of the 
j — > extrapolation and a better understanding of the 
volume dependence and appropriate normalisation of our 
results. The BEC window would be expected to reappear 
for smaller TO^/m p . 

While we have clearly defined the finite-temperature 
deconfinement transition at fi = 0, and find clear evi- 
dence of a deconfinement temperature that decreases as 
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/i increases, the exact nature and location of this transi- 
tion at large /z remain elusive. In order to pin down this 
transition, and also to precisely locate the superfluid-to- 
normal transition, we need to perform fine temperature 
scans by varying N T at fixed chemical potential. This 
is currently underway. We are also studying the static 
quark potential, which should give further insight into 
the nature of this transition. 

For the first time in this paper we have attempted 
to calculated renormalised energy densities via an esti- 
mate of Karsch coefficients obtained from simulations on 
anisotropic lattices. We find the resulting corrections to 
our earlier results are substantial, and indeed strongly 
suggest the quark contribution e g (/x) is negative, imply- 
ing that the physical requirement e q + e g > arises from 
a cancellation between terms of opposite sign. Consider- 
ably greater accuracy will be required, therefore, before 
we can contemplate eg. using lattice results as input for 
the solution of the Tolman-Oppenheimer-Volkoff equa- 
tions used in modelling relativistic stars. 

Wc find that the quark number susceptibility \q is re ~ 
markably independent of the temperature up to T ~ 
lOOMcV, and stays close to its noninteracting value in 
the quarkyonic region. Most strikingly, it shows little 
if any sensitivity to the deconfinement transition, which 
occurs at different chemical potentials for our 4 temper- 
ature values. 

The observation that Xq is n °t a proxy for the Polyakov 
loop L in the regime of high quark number density sug- 
gests the following conjecture. In the quarkyonic region 
fi < A* < M<2 the bulk obscrvablcs p, n q , Xq are approxi- 
mately equal to the free-field values pss, usb and Xsb- 
This is indicative of weakly self-bound quark matter, ie. 
with Ef = /U ~ &f. The transition at \i = [id (which co- 
incides with deconfinement as signalled by L ^ only in 
the limit T — > 0) is to a more strongly self-bound regime 
with Ef < kp. Since the quarkyonic phase is confining, 
we interpret weak self-binding as the quarks interacting 
via binary short-ranged interactions. For fi > [id, the in- 
teraction is screened and may not be much longer-ranged, 
but in this case deconfined quarks may interact with 
several other quarks in the vicinity leading to stronger 
binding. These considerations are related to interactions 
within bulk quark matter, involving quarks with all en- 
ergies less than the Fermi energy, which hence are not 
sensitive to temperature T. 

By contrast, the observed temperature-sensitivity of 
the Polyakov loop L (see Fig. [6]) suggests that in this case 
the relevant physics is associated with degrees of freedom 
close to the Fermi surface, which are readily thermally 
excited. These, of course, are the same degrees of free- 
dom relevant for transport. Our results contrast with the 
findings of analytic and numerical studies of QCD-likc 
theories at weak coupling in small volumes of character- 
istic scale R <C Aq [^S HI] > an d a recent study of cold 
dense QCD with heavy quarks [3Cj , both of which show a 
coincidence in the rise of L and Xq- This suggests that a 
full description of deconfinement at high baryon density 



requires a thermodynamic limit and light, mobile degrees 
of freedom. 

The main shortcoming of this study is that it has been 
performed with a single, relatively coarse lattice spacing. 
Although, as observed in Ref. [|[ , the main results are in 
qualitative agreement with the earlier results Q obtained 
on a coarser lattice with a = 0.23fm, we also observe sig- 
nificant quantitative discrepancies, and substantial lat- 
tice artefacts for fia > 0.75. To get this under control 
it will be necessary to repeat our simulations on a finer 
lattice. Thanks to the extensive investigation of param- 
eter space reported in Sec. [TH wc are in a good position 
to carry this out, and these simulation are underway. 

The large quark mass is another source of systematic 
uncertainty. It is clear that QC2D must be treated as 
a separate theory and cannot be viewed as an approx- 
imation to QCD - indeed, the differences between the 
two theories become most stark in the chiral limit - and 
there is hence no need to attempt to match quark masses 
to those in the real world. Still, many analytical results 
have been obtained in or near the chiral limit. Also, 
as already mentioned, we would expect a BEC region 
to open up near [i for smaller values of m w /m p , and a 
fuller understanding of the BEC-BCS crossover would be 
valuable. For all these reasons, simulations with smaller 
quark masses would be of great interest, and such simu- 
lations are underway. 

In addition to the quantities considered here, wc are 
in the process of computing the Landau-gauge gluon and 
quark propagators. This will allow us to check the as- 
sumptions involved in model solutions of the supcrfluid 
or superconducting gap equation, and may form a di- 
rect link with functional methods such and the functional 
renormalisation group and Dyson-Schwinger equations. 
These do not suffer from the sign problem, but rely on as- 
sumptions regarding the form of propagators and higher 
order vertices. This will be addressed in a forthcoming 
publication. 
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